﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class glDRIVE
{
    /*
    函数 gl.gear
    积分刚性方程组的Gear方法
    ss 计算雅可比矩阵的函数名。
    f 计算方程组中各方程右端函数值的函数名。
    参数 a: 积分区间起始点。
    参数 b: 积分区间终点。
    参数 hmin: 积分过程中允许的最小步长。
    参数 hmax: 积分过程中允许的最大步长。
    参数 h: 积分的拟定步长。hmin<<h<hmax。
    参数 eps: 误差检验常数。
    参数 n: 方程个数，也是未知数个数。
    参数 y0: y0[n]n个未知函数在起始点处的函数值。
    参数 k: 拟定输出的积分点数。
    参数 t: t[k]返回实际有效输出点（包括起始点）的自变量值。
    参数 z: z[n][k]返回实际有效输出点处的未知函数值。
    返回值 
    函数返回实际输出的积分点数。在函数返回之前，会给出下列相应信息供参考:
    全区间积分成功。若此时输出点数不够，可增大积分区间终点值
    步长小于hmin，精度达不到，积分停止(前输出点有效)
    阶数已大于6，积分停止(前输出点有效)
    对于h>hmin校正迭代不收敛，积分停止(前输出点有效)
    精度要求太高，积分停止(前输出点有效)
    */

    public static unsafe string drive_gear()
    {
        int i, j, k, m;
        double a, b, hmax, h;
        double[] y = new double[3];
        double[] t = new double[30];
        double[,] z = new double[3, 30];
        double[] hmin = new double[4] { 0.0001, 0.0001, 0.00001, 0.00001 };
        double[] eps = new double[4] { 0.0001, 0.00001, 0.00001, 0.000001 };

        a = 0.0;
        b = 1.0;
        h = 0.01;
        hmax = 0.1;
        gl.ss_x_ya_n_da = gears;
        gl.f_x_ya_n_da = gearf;
        string rs = "";
        for (k = 0; k <= 3; k++)
        {
            y[0] = 1.0;
            y[1] = 0.0;
            y[2] = -1.0;
            m = gl.gear(a, b, hmin[k], hmax, h, eps[k], 3, y, 30, t, z);

            rs += gl.html_table("k,h,hmin,hmax,eps", new double[] { k, h, hmin[k], hmax, eps[k] });
            rs += gl.html_table("实际有效输出点（包括起始点）的自变量值 t", t);
            rs += gl.html_table("实际有效输出点处的未知函数值 z", z);
        }
        return rs;
    }

    // 计算方程组各方程右端函数值
    private static unsafe void gearf(double t, double* y, int n, double* d)
    {
        //t = t; n = n;
        d[0] = -21.0 * y[0] + 19.0 * y[1] - 20.0 * y[2];
        d[1] = 19.0 * y[0] - 21.0 * y[1] + 20.0 * y[2];
        d[2] = 40.0 * y[0] - 40.0 * y[1] - 40.0 * y[2];
        return;
    }

    // 计算雅可比矩阵
    private static unsafe void gears(double t, double* y, int n, double* p)
    {
        //t = t; y[0] = y[0];
        p[0 * n + 0] = -21.0; p[0 * n + 1] = 19.0; p[0 * n + 2] = -20.0;
        p[1 * n + 0] = 19.0; p[1 * n + 1] = -21.0; p[1 * n + 2] = 20.0;
        p[2 * n + 0] = 40.0; p[2 * n + 1] = -40.0; p[2 * n + 2] = -40.0;
        return;
    }
}